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Abstract. We present an elementary and self-contained account of the 
analogies existing between classical diffusion and the imaginary-time evo- 
lution of quantum systems. These analogies are used to develop a new quan- 
tum simulation method which allows to calculate the ground-state expecta- 
tion values of local observables without any mixed estimates nor population- 
control bias, as well as static and dynamic (in imaginary time) response 
functions. This method, which we name Reptation Quantum Monte Carlo, 
is demonstrated with a few case applications to "^He, including the calcu- 
lation of total and potential energies, static and imaginary-time dependent 
density response functions, and low-lying excitation energies. Finally, we 
discuss the relations of our technique with other simulation schemes. 



1. Introduction 

The theory of stochastic processes plays an important role in modern devel- 
opments of quantum mechanics both as a deep — and possibly not yet fully 
exploited — conceptual method [1, 2] and as a powerful practical tool for the 
computer simulation of interacting quantum systems [3]. Quantum Monte 
Carlo simulations mainly rely on the static properties of one kind or the 
other of random walk that is used to sample the ground-state wave-function 
or finite-temperature density matrix of a system. Comparatively minor at- 
tention has been paid so far to the dynamical properties of the random 
walks used in quantum simulations. The main interest in these properties 
stems from the study of excitation energies, a notoriously difficult and ill- 
conditioned problem. These dynamical properties, also determine the mag- 
nitude of autocorrelation times which are necessary to estimate statistical 
errors. 



2 



STEFANO BARONI AND SAVERIO MORONI 



The purpose of these lectures is to provide an elementary and self- 
contained presentation of the deep relations existing between diffusion in 
classical systems and the imaginary-time evolution of quantum systems, 
and to develop some ideas based on these relations [4] which lead to a new, 
promising technique for performing quantum simulations. For reasons which 
will become apparent in the following, we name this technique Reptation 
Quantum Monte Carlo (RQMC). The main features of RQMC are that 
it is based on a purely diffusive process without branching, that ground- 
state expectation values of local observables can be evaluated without any 
mixed estimates, and that static and dynamic (in imaginary time) response 
functions are natural by-products of the ground-state simulation. Although 
some of the ideas presented here are not new [4] , our method does not suffer 
from the drawbacks which affected previous their implementations. 

2. From classical diffusion to quantum mechanics 

The simplest phcnomenological description of classical diffusion is given by 
the Langevin equation: 

dx = f{x)dT + d^, (1) 

where x is some generalized coordinate describing our system, f{x) = 
— ^g^^ is the force acting on it — which we suppose to be derivable from a 
potential — r is time and ^(r) is a Wiener process: {d^) = 0; ((d^)^) = 2dT. 
Although the random walk generated by the Langevin equation, Eq. (1), 
can be given a mathematically unambiguous meaning in the continuous 
limit, it is simpler — and more useful in view of applications to quantum 
simulations — to specialize to a given discretization of time: Tk = k x e. The 
corresponding random walk is then described by the discrete Markov chain: 

Xk+i=Xk + ef{xk)+^k, (2) 

where the ^'s are uncorrelated Gaussian random variables of zero mean, 
i^i) = 0; and variance 2e, (CiCj) = 2e5ij. 

2.1. THE FOKKER-PLANK EQUATION 

The time evolution of the probability distribution for the variable x is 
given by the master equation which, in the one-dimensional case, reads 
(the generalization to many dimensions is straightforward): 

P(x, r + e) = J Weix, y)P{y, r) dy, (3) 

where 
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is the conditional probability that the system is in configuration x at time 
r + e, given that it is found in configuration y at time r. When the con- 
ditional probability is independent of r (as it is in the present case) the 
corresponding Markov process is said to be homogeneous. In order to con- 
vert the master equation, Eq. (3), into a differential equation, we perform 
a Taylor expansion of its right-hand side in powers of the time step e: 

P{x, T + e) = [5{x-y- ef{y) - O^-^Piv, r) dy. (5) 

V47re J 

Taking into account that ^ ^ -^e, we now formally Taylor-expand the S 
function in powers of (ef{y) + 0> obtain: 

P(x,r + 6) = |p(y,r)(^X^^<5W(x-y)x((6/(y)+0")) dy, (6) 

where (•) indicates the Gaussian average of the polynomials in ^, and 5^'^^ 
is the n-th derivative of the 5 function: 

{{ef{y) + 0") ^ ^ r (efiy) + O^e'^d^ 
V47re J-<x> 

= , (7) 

Hn being the Hermite polynomial of order n [5]. To linear order in e, the 

first few values of these Gaussian integrals are: {ef{y)+(,) = ef{y), {{ef{y) + 
^)2) = 2e + (e^), {{ef{y)+(,f) = O (e^). Integrals of the derivatives of the 
6 function are given by: / g(y)6^"\x — y)dy = g^'"'\x). Using this relation, 
Eq. (6) can be recast as: 

P{x,r+e) = P{x,r)+e (^-^ (/(x)P(x,r)) + ^-^^^ +0 (e^) . (8) 

Eq. (8) is the discrete-time version of the Fokker-Planck [6] equation which, 
in the continuous limit, reads: 



dP{x,T) _ d^P{x,T) _ d 
dr dx'^ dx 



-{f{x)P{x,T)). (9) 



The fact that this equation is first order in time is strictly related to the 
Markovian character of the random walk, Eq. (2). It is easily verified that 
Ps{x) oc e"^*^^) is a stationary solution of the Fokker-Planck equation, Eq. 
(9). Using Eq. (8) one sees that the stationary distribution of the discrete 
random walk, Eq. (2), differs from Pg by terms of order e. 
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The conditional probability, We, defined in Eq. (4) can be seen as the 
exact probability that the system described by the discrete Markov chain, 
Eq. (2), makes a transition from configuration y to confignration x during 
the time step e; alternatively, it can be seen as an approximation to the 
transition probability in the continuous limit, correct to order O (e^) . In the 
continuous limit the exact conditional probability is defined by the relation: 

P{x, t) = j W{x, y- T)P{y, 0)dy. (10) 

By inserting this definition into Eq. (9), one sees that W{x, y; r) satisfies Eq. 
(9) with respect to x, subject to the boundary condition: W{x, y;0) = S{x — 
y), i.e. yV{x,y;T) is the Green's function of the Fokker-Planck equation, 
Eq. (9). The Markovian character of random walk, Eq. (2), allows one 
to define a simple functional-integral representation for W{x, y; r) which 
closely resembles the path-integral representation of the Green's function of 
the time-dependent Schrodinger equation [7] . As we will see in the following, 
this resemblance is by no means accidental nor superficial. 

Let X]\f = {xo,xi, ■ ■ ■ ,xi\f} be a given random walk generated by Eq. 
(2). Because of the Markovian character of the chain, the corresponding 
probability density, "PefXjv] satisfies the relation: 

Ve[XN] = Proh[xiNe) = XN;x{iN-l)e)=XN-i;---;x{0)=xo] 
= We{xN,XN-i)Proh[x{{N - l)e) = xn-V,- ■ ■ ■,x{0) = xq] 

= We{xN,XN-l)Ve[XN-l], (11) 

where x{t) is the configuration of the system at time r. By iterating this 
equation N times, one obtains: 

Ve[XN]=mixN,XN-l)m{xN-l,XN-2)---m{xi,Xo)Pixo,0). (12) 

The probability density that the system is in configuration xn at time 
r = Ne is obtained from VeiXj^] by integrating out the N variables, 

{xq, ■ ■ ■ ,XN-l}- 

P{xN,Ne) = J Te[{xo,--- ,XN}]dxo---dxN-i (13) 

= j yVe{xN,XN~l) ■ ■ ■yVe{xi,Xo)P{xo,0)dxo - ■ ■ dXN-l- 

By comparing Eq. (10) with Eq. (13), one obtains the desired functional- 
integral representation for W: 

W{x,y;T) = 

J yVe{x,XN-l)y\^e{xN-l,XN-2) ■ ■ ■y\^e{xi,y) dxi---dxN-l, (14) 
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where e = t/N . The above representation is exact for the discrete Markov 
chain described by Eq. (2). In the continuous hmit (e — > 0) it is understood 
that it holds by letting iV = r/e — oo, while keeping r fixed. 

2.2. THE CLASSICAL-QUANTUM MAPPING 

In order to establish the link between classical diffusion and imaginary- 
time quantum evolution, we define a wave-function, $(x,r), through the 
relation: 

P{x,t) = ^q{x)^{x,t), (15) 

where ^{){x) = ^J Ps{x) oc e""*^^)/^. By inserting Eq. (15) into Eq. (8) and 
dividing by ^q{x), we obtain: 

^{x, r + e) = (1 - eU) $(a;, t) + O (e^) , (16) 

where: 



|2 



^ = -^+V(x), (17) 



and 

V{x) 



1 (dv\^ _ 1 d'^v 
4 \dx) ~2dx^ 

1 #*o(x) (jg, 



^o{x) dx^ 

The continuous-time limit of Eq. (16) is formally equivalent to an imaginary- 
time Schrodinger equation, 

'-^^-nM.,r). (19) 

which could as well have been derived directly by inserting Eq. (15) into 
Eq. (9). 

The wave-function $o is a solution of Eq. (19), which means that it 
is an eigenfunction of the time-independent Schrodinger equation, corre- 
sponding to a zero eigenvalue. A general theorem of quantum mechanics 
states that the ground-state eigenfunction of a Schrodinger equation with 
a local potential is non-degenerate and node- less [8]. Orthogonality with 
respect to excited-state wave-functions implies that the ground state is the 
only node-less eigenfunction. We thus arrive at the conclusion that all the 
excited-state energies are strictly positive and that ^o{x) is the only time- 
independent solution of Eq. (19) or, equivalently, that Ps{x) oc e"''^^^ is 
the only stationary solution of the Fokker-Planck equation (9) [9] . Further- 
more, any solution of Eq. (9) would tend to Ps for large times, irrespective 
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of the initial condition. If the spectrum of 7i has a gap, the approach to 
equihbrium is then exponentially fast. This is easily seen by simple inspec- 
tion. Consider a system whose probability distribution at time r = is 
0) and its expression in terms of the eigenfunctions and eigenvalues of 
H, which we indicate by and £n- 

P{x,Q) = ^Q{x)Y,Cn^n{x). (20) 
n 

The time evolution of P is readily derived from the time evolution of the 
$'s: 

P(x, r) = c^^aixf + c„e-^"^$o(a;)$n(x)- (21) 

The normalization of P{x,G) and the orthonormality of the $'s imply that 
Co = 1 and that the norm of P{x, r) is conserved. The above equation 
shows that the thermalization time — i.e. the time necessary for the system 
to reach equilibrium — is tqk, 

The fact that Pg is the only stationary solution of the Fokker-Planck 
equation, Eq. (9), implies that classical expectation values over Pg — or, 
equivalently, quantum expectation values over ^o{x) — can be expressed 
as time averages over the random walk, x{t), generated by the Langevin 
equation, Eq. (1): 

/ Ps{x)A{x) dx= {^o\A\^q) = lim lim - / A{x(t)) dr. (22) 

J T-^oo t-^oo t Jt 

A comparison between Eq. (3) and Eq. (16) allows one to establish a 
relation between the transition probability of the classical random walk. 
We, and the propagator of the quantum system associated with it: 

We{x,y) = Mx)g{x,y;e)/My) + 0(e^) , (23) 

where: 

g{x,y;T) = {x\e-^^\y). (24) 

The fact that the error in Eq. (23) is indeed of second order in e, and not 
higher, can be proved by pushing the Taylor expansion of Eq. (6) to second 
order in e and by noting, for instance, that the second-order term would 
give rise to a non-hermitian contribution to Q. Using Eqs. (12) and (23), 
the probability density for the random walk X^r can be easily expressed in 
terms of a product of quantum propagators, Q. Assuming that the system 
is at equilibrium at r = — i.e. that the probability distribution for xq is 
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Psixo) — the probability distribution for the random walk is: 

Ve[XN] = ^o{xN)Q{xN,XN-i;e)---Gixi,xo;€)^o{xo)+0{e) 

V[Xn] 

= V[X]x Q.IXn], (25) 

where Qel-'^iv] = 1 + O (e) is a time-discretization correction factor. The 
Markovian character of the random walk implies a simple composition law 
for the V probability distributions: 

P[{xo, • • • , Xjfc, • • • , Xn}] = V[{Xo, ■ ■, Xk}] X 

V[{xk,---,XN}]/Ps{xk). (26) 

Inspection of Eq. (25) shows that the probability distribution 7^[-X'] is in- 
variant under time reversal: 

V[X]=V[X], (27) 
where X = {xn, ■ ■ ■ , xq} is the path obtained from X under time reversal. 

2.3. CLASSICAL VS. QUANTUM TIME CORRELATION FUNCTIONS 

Eq. (22) expresses the identity between classical thermal expectation val- 
ues of local observables and quantum expectation values of the same ob- 
servables calculated over the ground state of a suitably defined auxiliary 
system. Furthermore, these expectation values can be expressed in terms 
of time averages over a random walk. In the following we will see how this 
result can be extended to time correlation functions. It will result that 
(real) time correlation functions calculated over the random walk coincide 
with the ground-state imaginary-time correlation functions of the auxiliary 
quantum system. 

The time auto-correlation function of an observable A{x) is defined as: 

{A{t)A{0)) = j A{x)A{y) Prob[a;(0) = y\x{T) = x] dxdy. (28) 

If the stochastic process is stationary [i.e. if P{x,t) = Ps{x)), then one 
has: 

Prob[a;(0) = y; x(r) = x] = yV(x, y; t)Ps {y) 

= $o(ic)g(x,y;r)$o(y)- (29) 
By inserting this relation into Eq. (28) , the auto-correlation function reads: 

{A{t)A{0)) = J A{x)A{y)Mx)G{x,y;T)^o{y)dxdy 

= ($o|e^Me-^M|$o> 

= (^>o|^(-ir)^(0)|$o), (30) 
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where we have used the fact that e^'^|$o) = |^o) and the definition of a 

time-dependent operator in the Heisenberg representation: A{t) = e^^^A 

Q—iHt 

Eq. (30) can be used to express the autocorrelation time of the ob- 
servable A in terms of the spectral properties of the auxiliary quantum 
Hamiltonian, Ti. The autocorrelation time is defined as the time integral of 
the normalized autocorrelation function: 

'^ = Jo {A-) - {A}- ^^'^ 

^ V ' 

The integrand in Eq. (31) is defined in such a way that c^(0) = 1 and 

c^(oo) = 0- When c_4 is exponential, is simply its decay constant: 
ca{t) = e""^/"^-^. In practical simulations, r4 determines the statistical er- 
rors in the measure of A: 

M+N 



E A{xi)±JAA^^ (32) 

i=M+l * 

where M is chosen large enough so that equilibrium is reached (i.e. M > 
^), and AA^ « jv^Z)(-4(xj) — {A})'^. In order to proceed further, we 
first consider the spectral resolution of the quantum propagator, Q: 

g{x, y;T) = Y^ e-^""$n(x)$„(y). (33) 

n 

By inserting Eq. (33) into Eq. (30), the auto-correlation time can be easily 
cast into the form: 

1 K^^4^ 
^■A = ^ 2^ g • (34) 



3. From quantum mechanics back to classical diffusion 

The purpose of many quantum simulation techniques is to study the ground- 
state properties of a system whose Hamiltonian is 

and whose (unknown) ground-state wave-function and energy we indicate 
by \I'o and Eq, respectively. In the variational Monte Carlo method (VMC), 
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an approximate wave-function, is used to generate a random walk ac- 
cording to the discrete Langevin equation, Eq. (2), with 

m^fvMci^) = -^(-iog$o(x)2) 

~ $o(x) dx ' ^'^^^ 

and the ground-state expectation value of the operator A is estimated 
through Eq. (32), where A{x) = ^^^A^o{^)- Systematic errors in Eq. 
(32) depend on the discretization of time and are of order e. They can 
be eliminated in principle by enforcing the detailed-balance condition on 
the master equation, Eq. (3) [10]. A variational upper bound, Sq, to the 
ground-state energy can be estimated from Eq. (32): 



where 



= ^7-Ti^^o(x) (38) 

is the so-called local energy. In the following we will show how the dynamical 
properties of the random walk (1) can be used to systematically improve 
upon this VMC procedure and to estimate, exactly within statistical noise, 
the ground-state properties of quantum systems. 

Let us first observe that the trial wave- function, <I>o, implicitly defines a 
reference (unperturbed) Hamiltonian, Ho, whose exact ground state is $o- 
The potential function which defines Hq is simply obtained by inverting 
the time-independent Schrodinger equation: 

A comparison between Eq. (39) and Eqs. (17,18) shows that the potential 
obtained from the inversion of the Schrodinger equation coincides with the 
effective quantum potential resulting from the classical-quantum mapping 
discussed in Sec. 2.2, provided that the classical random walk, Eqs. (1,2) 
is driven by the VMC force defined in Eq. (36). The original Hamiltonian, 
Eq. (35), can then be cast into the form: 

H = -^ + Vo{x) + £{x)-So. (40) 
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By construction, $o is the ground-state of TL, and £{x) = Eq if $o is the 
ground state of H. When $o is not the ground state of iJ, the local energy 
£{x) is not a constant, and the ground-state properties of H can in principle 
be obtained by perturbation theory with respect to A7i. 

Let us calculate the corrections to the ground-state energy, £q, up to 
second order in A7i. The first-order correction vanishes by construction, 
given that ($o|^"|^o) = £o- The second-order correction is given by: 

A£:f = -i:i»^. (41) 

By comparing this expression with Eq. (34), the second-order correction to 
the ground-state energy can be expressed in terms of the fluctuations of 
the local energy and of its auto-correlation time: 

a4') = Te^£\ (42) 

In the following, we show how this relation between perturbative corrections 
to the ground-state energy and (integrals of) auto-correlation functions of 
the local energy along the random walk can be generalized to arbitrary or- 
der. The resulting perturbation series can be effectively summed to infinite 
order by an appropriate re-sampling of the random walks generated by the 
VMC Langevin equation, Eqs. (1,2,36). 



3.1. STOCHASTIC PERTURBATION THEORY 

Given the Hamiltonian H, Eq. (35), and a trial wave-function, which 
we suppose to be non-orthogonal to its ground state, ^'o, the ground-state 
energy of H can be expressed as: 

^0= lim = - lim #log($o|e-^1$o) (43) 

= - lim ilog($o|e"-^^|^>o). (44) 

r— >oo 7" 

These equations are easily demonstrated by expressing $o as a linear com- 
bination of the eigenfunctions of H, {^'n}: = J2n Cn^n-, and by inserting 
this expression into Eqs. (43,44): 

log(cl>o|e-^^|<l>o) = -E^T + log (4) + O (e-(^i-^°)^) . (45) 
The following basic identity holds: 

Zo = ($o|e-^n$o) = y e-«^[^lp[X] 

- (e-^[^l)> (46) 



REPTATION QUANTUM MONTE CARLO 



11 



where 7-*[X] is the probabihty distribution for the random walk X, Eq. (25), 
"DlX] = dxo ■ ■ - dxN, and the action S[X] is a functional of the random walk 
which in the continuum limit coincides with the time integral of the local 
energy: 

TV 

S[X] = 6^^(x,) + 0(e) 

i=l 

« r £{x{T'))dT' . (47) 
io 

Eq. (46) is a generalization of the Feynman-Kac formula [11] and can be 
demonstrated as follows. 

Zo = j $o(a;)G(x,y;r)#o(y) dxdy, (48) 

where 

G(x,y;r) = (x|e-^ny) (49) 

is the imaginary-time propagator of the full Hamiltonian. We now break 
the propagator in Eq. (48) into the product of N short-time propagators: 

2'0 = j po(xo)G'(xo, xi; e) • • • G{xn-i,xn] e)^o(a;jv) X>[X]. (50) 

Eqs. (46) and (50) can be seen as a definition of the action: 

V[X]=V[X]e-^^^l (51) 

By using the Trotter formula, 

G{x, y- e) = g{x, y; e)e-^ ^(^) + O {f) , (52) 

one sees that the action defined by Eq. (51) can indeed be expressed by the 
time integral given by Eq. (47). 

Zq plays the role of a pseudo partition function, in the sense that the 
energy of the system — as well as other observables, as we will see — can be 
calculated by differentiating it much in the same way as one would do in 
classical statistical mechanics. In particular, the two expression (43) and 
(44) for the ground-state energy correspond to the zero-temperature limits 
of the internal and free energies, respectively. By inserting Eq. (46) into Eq. 
(44), one could derive a systematic perturbative expansion of the ground- 
state energy in powers of £. Each term of the series is basically a cumulant 
of the action which in turn can be seen as the integral of a suitably defined 
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connected time correlation function of the local energy, calculated along the 
Langevin random walk. This kind of stochastic perturbation theory can be 
effectively carried on to infinite order by inserting Eq. (46) into Eq. (43). 
The final result reads: 

(£:(a;(r))c-'5W 
En = lim ^ 



^ ^lim {(f(x(r)))), (53) 

where the double bracket ((•)) indicates that the average over the random 
walks {quantum paths) is re-weighted by the exponential of the action, Eq. 
(47). 

Eq. (53) can be turned into an algorithm for calculating the ground-state 

energy. The calculation would proceed as in a standard VMC simulation, 

with the difference that the local energy must be weighted with the expo- 

— ^ 

nential of the action, e"*-^ , calculated along a segment of the random walk 
long T and ending at the time when the measure is taken. This algorithm, 
which was first proposed in Ref . [4] , is bound to fail in all those case where 
the number of particles is so large or the quality of the trial wave-function 
is so poor that the fluctuations of the action make the weighting procedure 
impractical. As a matter of fact, the pure- diffusion Monte Carlo (PDMC) 
of Ref. [4] has never been applied but to very simple quantum systems. 
In the next session we will show how the fluctuations of the action can be 
effectively dealt with through a new algorithm based on importance sam- 
pling. Before doing this, we want to show how the ideas exposed so far 
can be exploited to calculate general ground-state expectation values and 
response functions. 

3.2. CALCULATION OF OBSERVABLES 

Other observables can be calculated along similar lines starting from the 
Hellman-Feynman theorem [12]: 



(54) 

A=0 



dX 

where Ex is the ground-state energy of a system whose Hamiltonian is 

Hx = H + XA. (55) 
By using Eq. (44), Ex can be put in the form: 

= - jin^ z log (^"^^ ) ' (56) 



T— >00 f 
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where 

/o 



Sx[X] = r{S{x{T')) + \A{x{t'))) dr'. (57) 
Jo 

The expectation value of A can then be expressed as: 

l/oM(x(r'))dr'e-/o^^(-')<^r'^ 



r JO 

(*o|^|*o) = lim 



3.3. RESPONSE FUNCTIONS 

A simple extension of the ideas used in the previous section leads to a 
technique for evaluating response functions. Let us suppose that the system 
is coupled to a set of external variables, {Aj} through the operators {Ai}: 

H{x} = H + Y,)^iAi. (59) 

i 

In the case of an external potential, for instance, the index i labels the 
coordinate, x, Aj is the potential itself, V(,xt{x)-, and Ai is the particle- 
density operator, n{x). Wc define a generalized susceptibility, Xi ^-s the 
derivative of the expectation value of one of the A operators with respect 
to one of the A's: 

d{A) 

Xij - 



dXidXj 



(60) 



By using Eq. (56), x can be put in the form: 

' r Ai{T')dT' r Aj{T')dT' 

\Jo Jo 



Xij = - lim - 

T—^00 T 







2' 

AiT')dT' 



= - Hm ^ ((^J^ dn £ dr2iAiiTi)Aj{T2) - Ai ^,-))) , (61) 

where A is the time average of A over the random walk: A = ^ Jq A{T')dT'. 
We now split the domain of integration [0 < ri < r; < T2 < r] into two 
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sub-domains [ri < T2] and [t2 < ti], and change the variables of integration 
{ti,T2} — {ti,T2 — Ti} and {ti,T2} — ^ {n — r2, T2} in the two sub-domains 
respectively. In the large r limit, the susceptibility can then be cast into 
the form: 



The symmetrized time correlation function Cjj + Cji is a functional of the 
path, whose time integral (i.e. whose a; = Fourier component) provides 
an estimator for the static response function. This argument can be easily 
generalized to show that other (uj 7^ 0) Fourier components of the same 
time correlation function provide suitable estimators for the dynamic sus- 
ceptibility. 

4. The algorithm 

In the previous section we have shown that the calculation of ground-state 
expectation values and response functions can be reduced to the weighted 
average of suitable estimators over the space of random walks. A straight- 
forward evaluation of these averages over a Langevin trajectory would be 
very inefficient because the weight, e , may vary a lot. This problem 
can be solved by converting the weighting of the quantum paths into a 
re-sampling of them, through an appropriate Metropolis algorithm [13]. As 
a by-product, the use of the Metropolis algorithm allows one to reduce 
the systematic errors due to the discretization of time in an efficient and 
convenient way. 

Let G^'"^ be any approximation to the full imaginary-time propagator, 
Eq. (49), correct to order n in e, and P(")[X] the corresponding approx- 
imation to P[^], Eq. (50), which will be affected by errors of order n. 
Let us also define the corresponding approximations for the unperturbed 
propagator, Eq. (24), path probability distribution, P^"), and time- 

discretization correction factor, Q{^\ Eq. (25), and action, S^'^^[X], Eq. 
(51). For instance, one can take: 




(62) 



where the time auto-correlation function is defined by: 




g(2)(x,y;e) 

p(2)[X] 



^Q{xt^)Q^ '{xN,XN-i;e)---Q^ ^(xi,a;o;e)$o(a;o) 
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= V[X]+o(€^), (64) 



5(2) [X] = e-t Er=i(£(-*)+^(-*-i)) = S[X] + O (e2) . 

Any ground-state expectation value or response function can be put in the 
form: 

/p(^)mF|xi^ 

\\ L J// /p(«)[X] T>[X] ^ ' ^ ^ 

In order to calculate the above expectation value through the Metropolis 
method, it is necessary to construct a Markov chain of random walks so 
that the corresponding transition probability, W[Y X], satisfies detailed 
balance: 

W[Y ^X]x P(") [X] = W[X ^Y]x P(") [Y] (66) 

In the Metropolis algorithm, the transition probability W is split into the 
product of an a-priori transition probability, W*^, times a probability. A, 
that a move proposed according to W'' is accepted: W[y ^ X] = W*^[y <— 
X] X A[Y <— X]. Detailed balance, Eq. (66), can be satisfied by choosing 
[14]: 

A[Y ^X]= min(l, R[Y ^ X]), (67) 

where 

w°[x^y]xp(")[y] 
^ - wo[Y^x]xp(n)[x] ■ ^^^^ 

Given a quantum path, X = {xq, ■ ■ ■ , xn}, we propose a new path, Y, by 
propagating the random walk forward by M < A?^ steps, according to the 
VMC Langevin equation, Eqs. (2,36): Y = {xm, ■ ■ ■ ,xn, ■ ■ ■ ,xm+n}- The 
corresponding a-priori transition probability is: 

W^[Y^X] = W,{xN+M,XN+M-l)---We{xN+l,XN) (69) 
= V^'^^[{XN, ■ ■ ■ ,XN+M}]Qi''^[{xN, ■ ■ ■ ,Xn+m}]/Ps{xn), 

where we have used Eq. (25). By inserting this expression for the a-priori 
probability into Eq. (68), the latter can be cast into the form: 

^.y ^X]= [{^M, • • • , xq}] X [{xm, ■ ■ Xo}]/PsixM) 

^ V^''^[{xn,---,Xn+m}]x Qi^^[{xN,---,XN+M}]/Ps{xN) 



X 



e-«5^"^W 

pH[{xo, • • ■,xm}] X — (70) 



With the aid of the composition law for the random-walk probability dis- 
tributions, Eq. (26), and of the time-reversal property, Eq. (27), the above 
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equation finally reads: 

K[Y ^X]= e-(5^"'[ri-5'">[x]) ^ Q^W, • • • , xp}] ^^^^ 

Qr [{^N, - ■ ■ ,XN+m}] 

Notice that the ratio of the Q's goes to one in the continuous-time limit 

(e 0) and that for finite e it can be explicitly calculated, providing thus 
a systematic way for reducing to any desired order the time-discretization 
errors. 

The dynamical variables of our simulations are quantum paths (or seg- 
ments of random walks) which can be formally associated with polymers, 
much in the same spirit as this is done in path-integral simulations. The 
polymer dynamics which would correspond to our algorithm is known in 
the literature as reptation [15]. For this reason, we name our algorithm 
Reptation Quantum Monte Carlo (RQMC). 

The practical implementation of RQMC is extremely simple, at the level 
of a VMC simulation. The algorithm can be summarized as follows: 

1. Using Eqs. (2) and (36), generate a random walk long enough that its 
end point, x{T), is distributed according to $q. 

2. Generate a further segment of random walk corresponding to the time 
interval [T,T + t]. t = Ne should be large enough that the limit 
r — ^ oo in Eqs. (53), (58), and (61) is reached to the desired accuracy. 
Set X = {xo, xi, • • • , xn} ^ {x{T),x{T + e), • • • , x{T + r)}. 

3. Select a 'direction of time' (forward or backward) with equal probabil- 
ity. If the choice is backward, set X X = {xn,xn-i, • • • ,xo}. This 
step is necessary to ensure the micro-reversibility of the algorithm. 

4. Generate a segment of random walk corresponding to the time interval 
[T+T, T+T+S] and set Y = {xm, xm+i, • • • , xn+m} ^ {x{6),x{T+5+ 
e), - ■ ■ ,x{T+T+S)}. The value of 6 is sampled from an uniform deviate 
in the interval [0, A] whose width is chosen so as to minimize the auto- 
correlation times of the measured quantities. Sampling 5 instead of 
keeping it constant helps to avoid that the path remains occasionally 
stuck at a fixed position for a long time. 

5. Evaluate Ti\V <— X] according to Eq. (71). 

6. If R > 1, set X ^ y. If R < 1, set X ^ y with probability R. 

7. Accumulate the ground-state energy and other observables using ap- 
propriate estimators (see the next section for an optimal choice of the 
estimators). 

8. Go to 3. 

This algorithm has been preliminary tested for the hydrogen atom us- 
ing an approximate trial function. Exact results for the average of several 
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moments of electron-nucleus distance have been reproduced within a sta- 
tistical error pushed down to a small fraction of the difference between the 
exact value and the extrapolated estimate (i.e. twice the mixed average 
minus the variational average [16]). 

5. A case study of ^He 

We now discuss the calculation of several properties of superfluid ^He, with 
the purpose of showing that the method can be successfully applied to 
systems of actual physical interest. Based on the limited experience gained 
in this case study, we also present some performance comparisons with 
related techniques. 

5.1. DETAILS OF THE SIMULATION 

We consider Np = 64 ^Hc atoms interacting through a pair potential, as ob- 
tained from first-principles calculations [17]. The simulation was performed 
in a cubic box with periodic boundary conditions at the experimental equi- 
librium density, p = 0.02186 A^^. The trial function <I>o includes pair and 
nearly optimal three-body correlations [18]. This is a relatively good trial 
function. The variational energy is less than 0.3 K above the exact ground 
state energy, whereas the variational bias is larger than 1.1 K using pair 
correlations only. The quantities we compute are total and potential ener- 
gies, the imaginary-time correlations of the density fluctuation operators, 
Pq = Eiexp(-iq-ri): 

F{q,T) = {p,{T)plg{0))/Np, (72) 

and the diffusion coefficient of the center of mass motion, 

D{t) = ([rcM(r) - rcM(0)]2)7Vp/(6r). (73) 

The parameters of our simulations are as follows. The time step is e = 
0.001 K""^, which gives a systematic bias of the order of 10~^ K on the 
total energy. For the calculation of total and potential energy the length 
of the path was r = 0.4 K~^, corresponding to N = 400 time slices. The 
calculation of imaginary-time correlations over a significant range required 
the use of longer paths, up to A'^ = 700. The value of the energy resulting 
from the simulation with such longer paths confirmed that the finite-r bias 
in the results obtained with N = 400 was smaller than statistical errors. 
Note that the length of the path adversely affects the efficiency, because the 
relaxation time of the polymer in the reptation algorithm is proportional 
to N'^ [19]. The number of time slices of each reptation move is uniformly 
sampled between and 20, yielding an acceptance ratio of « 80%. 
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Figure 1. Average of the potential energy on individual time slices along the path. The 
statistical error on the central slices is ~ 0.03 K. This result was obtained using a trial 
function with pair correlations only: note that V converges to the same value given in 
Table 1, obtained using a trial function with pair and triplet correlations. 



Fluctuations in the average of the total energy are reduced using a sym- 
metrized form of Eq. (53), i.e. accumulating [£{x{t)) + £{x{Q))\ /2. Expec- 
tation values of local observables are computed averaging time integrals 
along the path, Eq. (58) and Eq. (61). Although these expressions are ex- 
act in the r — > oo limit, the contributions coming from the extrema of the 
path are clearly biased. For instance, the average of A in the initial or final 
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TABLE 1. Ground state energy, Eo, 
and potential energy, V, as computed 
from RQMC and traditional difltusion 
Monte Carlo runs of 3x10'^ Monte 
Carlo steps. Units are K. The length 
of the path in the RQMC calculation is 
T — 0.4 K^^, and the length of the for- 
ward walk for V in the diffusion Monte 
Carlo calculation is 0.2 K"'^. 

Eo V 

RQMC -7.4066(27) -21.644(15) 
BDMC -7.3902(15) -21.674(21) 



time slices {{A{x{0))) and {A{x{t)))) yields the mixed estimate: 

($o|^|*o)' ^ ^ 

On the other hand any time slice a distance r apart from the extrema, 
such that exp(— iff) |<I>o) — |^o)) gives an unbiased contribution to the 
time integral of A- Therefore it is convenient to restrict the time integral 
of A in Eqs. (58) and (61), to the inner section of the path, where the bias 
is reduced. For the potential energy and the imaginary-time correlations 
we exclude the contributions from 150 time slices on each side of the path. 
This is a rough estimate of the time it takes for the average potential 
energy to converge within a few htindrcdths K from its vahic at slice or 
A'',- (which is by Eq. (74) the mixed estimate) to the unbiased estimate. This 
is demonstrated by the average of the potential energy on individual time 
slices, shown in Fig. 1. Wc have not studied the corresponding convergence 
times for the density-density correlations. 

5.2. RESULTS 

Our results for the total and potential energies are listed in Table 1. The 
inter-particle potential adopted [17] overestimates the experimental bind- 
ing energy of —7.17 K because of the neglect of three-body forces (mostly 
triple-dipole repulsion). Also reported in Table 1 are the corresponding 
data obtained from a standard Branching Diffusion Monte Carlo (BDMC) 
calculation using the same time step and trial function. The small differ- 
ences between the results of the two algorithms have to be attributed to 
different time step biases. 
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Figure 2. Imaginary-time correlations of the density fluctuation operator. Averages are 
taken on the inner part of a path of length 0.7 K^^. The statistical error ranges from 
approximately 0.5% at t = up to 5% at r = 0.4 K~^. 



In figure 2 we show the density-density correlation function, F(q, r), for 
a few values of g, as obtained from a run of 10^ MC steps, which required 
about one week CPU time on a workstation. Note that at virtually no 
additional cost many more q vectors, belonging to the reciprocal lattice of 
the simulation box, could have been included in the calculation. 

F{q,T) is related to several quantities of physical interest [20], including 
the static structure factor, 



S{q)=F{q,0)= dwS{q,u), 
Jo 



(75) 
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Figure 3. Static structure factor at the five wave vectors listed in Fig. 2 (open circles). 
The solid line is the experimental S{q) measured by neutron scattering [21]. 

the static linear response function, 

poo roo 

Xiq) = -2 dr F{q, r) = -2 duj S{q, u:)/uj, (76) 
Jo Jo 

and the dynamical structure factor, 

F{q,T)= dw e-'^^ S{q,u;). (77) 
Jo 

Low moments of S{q,Lo) can be accurately extracted from the simu- 
lation data for F{q,T). The static structure factor, Fig. 3, and the static 
response. Fig. 4, compare very favorably with the experimental results, the 
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Figure 4- The static response function at the five wave vectors listed in Fig. 2 (open 
circles) . The solid line is the experimental result of Ref. [22] . 



discrepancy visible in S{q) at the smallest value of q being due to the finite 
temperature at which the measurements are performed. The /-sum rule, 
dF{q,T)/dT\r=o = J dio S{q,Lo)uj = q^, is also verified with high precision. 

Inferring dynamical properties from imaginary-time correlations, on the 
other hand, is much harder. Since the inverse Laplace transform (77) with 
incomplete and noisy data for F(q,T) is an ill-conditioned problem [23], a 
least approach to pin down a parametrized form for S{q, to) is doomed to 
failure. Additional constraints can be set on the solution S{q, u) by resorting 
to Maximum Entropy (ME) methods [23]. We follow the implementation 
used in Ref. [24] to process data for F{q, r) obtained at finite temperature 
with a PIMC simulation. The results are qualitatively similar to those ob- 
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Figure 5. Maximum Entropy reconstruction of the dynamical structure factor. In the 
inset the position of the maxima of the calculated S{q, uj) is compared with the experi- 
mental excitation spectrum [25]. 

taincd in Rcf. [24]. The ME reconstruction of S{q,uj), shown in Fig. 5, is 
too smooth and does not reproduce the sharp features exhibited by the 
experimental structure factor. Furthermore, the relatively poor quality of 
the available Monte Carlo data does not allow for a reliable estimate of the 
statistical uncertainty on the results [23]. Nevertheless we do recover some 
useful information on dynamical properties: the presence of a gap in the ex- 
citation spectrum is unambiguously revealed, and the position of the peak 
of the reconstructed dynamical response closely follows the experimental 
dispersion of the elementary excitations (see the inset of figure 5). 

We now outline the calculation of the superfluid density pg. Although 
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the value of ps/p is trivially one for pure bulk "^He in the ground state, 
interacting Bose systems in the presence of an external disordered potential 
undergo a zero temperature supcrfluid-insulator transition as the strength 
of the potential increases [26]. We thus consider a model system of static 
impurities in ^He. The external potential Vext is represented by attractive 
Gaussians, V"ea;t(r) = >lexp[— a(r — Rj)^], where the positions Rj of 
the impurities are placed randomly in the simulation box, A = — 50 K 
and a = 0.5 A~^. The trial function is multiplied by a one-body factor 
exp[— X^jj f{\ri — Rjl)] which tends to localize the ^He atoms around the 
impurities. No average over different realizations of disorder was performed. 

In a finite-temperature calculation with periodic boundary conditions, 
Ps can be estimated [27, 28] as ps/p = {w'^) /{GrNp), where the winding 



number w = J2i=i lo dr' '"^^T ' is the displacement of the center of mass 
of the system, and r is the inverse temperature. By taking the limit r — > oo, 
appropriate for a ground state calculation, the superfluid density at T = 
can be expressed in terms of the diffusion coefficient of the center of mass 
motion, ps/ p = lim^^co D(r), where D is defined in Eq. (73). 

The results reported in Fig. 6 show that the superfluid fraction, which 
is correctly one for the pure system, is indeed reduced in the presence of 
the impurities. Longer simulations would be needed to relate the depletion 
of the superfluid fraction induced by the disordered potential to changes of 
the excitation spectrum. 



5.3. COMPARISON WITH OTHER METHODS 



The RQMC method utilizes the Metropolis algorithm to sample an ex- 
plicitly known probability distribution, namely a discretized path integral 
expansion of exp(— ri?)$o which becomes exact in the limit t ^ oo and 
e ^ 0. Sampling a distribution, as opposed to carrying weights, avoids the 
difficulties associated with fluctuating weights which plague applications of 
'Pure Diffusion' [4] or 'Single Thread' [29] Monte Carlo. For instance, we 
were unable to get converged results for our 64 particle system with PDMC. 

The idea of sampling a path-integral representation of the imaginary- 
time evolution to compute ground state properties is not new [30]. For 
instance. Variational Path Integral [27] (VPI) uses the pair product ap- 
proximation to expand the many-body density matrix exp(— ri?) and the 
bisection method to sample the path, just like in the usual Path Integral 
Monte Carlo method [28]. 

RQMC features instead an expansion of the imaginary-time propaga- 
tor based on the Langevin dynamics generated by the trial function, and a 
reptation algorithm to sample the exponential of the resulting action. From 
the computational point of view, the advantage of this particular choice can 
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Figure 6. The diffusion coefficient of the center of mass motion. 



be understood in the limit of perfect importance sampling: if the trial func- 
tion is exact the local energy is a constant, and reptation moves consisting 
of an arbitrary number of time slices will be accepted with probability 1. 
Eventually, a very poor wave function (or equivalently a very large num- 
ber of particles) will force us to take extremely small reptation moves, and 
moves of the kind used in VPI will become more efficient [28] . VPl has been 
implemented [31] for the simulation of superfiuid ^He at T = 0. According 
to the author of the VPI calculation, for the system size considered here 
RQMC is considerably more efficient. 

Sampling an explicitly known distribution is to be contrasted to the 
standard BDMC, which samples an unknown distribution (i.e. the mixed 
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distribution ^'o^o) obtained from the asymptotic solution of a differential 
equation [32, 33]. BDMC is designed to compute efficiently the total energy; 
however it introduces a population-control bias, yields a mixed estimate for 
operators not commuting with H, and does not retain direct information on 
the imaginary-time correlations. This information can be retrieved through 
the 'forward walking' technique [33, 34, 35] and used to correct both the 
population control and the mixed estimate biases, but at the price of ad- 
ditional statistical fluctuations. We believe that RQMC will turn out to 
be advantageous over BDMC in those cases where recovering information 
from fluctuating weights becomes too noisy. 

The results for the total energy shown in Table 1 are obtained with 
the RQMC and the BDMC methods using the same time step, number of 
particles and trial function. Prom the estimated statistical error we infer 
that BDMC is roughly 3 times faster than RQMC for the calculation of 
the total energy. Also listed in Table 1 are the unbiased estimates for the 
potential energy. In this case the BDMC algorithm with forward walking 
(implemented in the "backward storing mode" described in Ref. [33]) turns 
out to be roughly two times slower than RQMC. Obviously, factors 2 or 3 
for a couple of observables in a particular physical system are not a conclu- 
sive assessment of the relative performance of two algorithms. However the 
resulting factor 6 in the relative improvement of RQMC when BDMC has 
to be complemented with forward walking suggests that, whenever explicit 
information on imaginary-time correlations is used, RQMC is likely to be 
competitive or better. 

6. Conclusions 

The most attractive feature of the RQMC method is that it uses the dy- 
namical properties of the random walk in a way that is directly related to 
the imaginary-time properties of the physical system. A previous imple- 
mentations of similar ideas (PDMC, [4]) was based on re- weighting and the 
scope of its applications was severely restricted by the fluctuations of the 
weights. Wc have shown that, by simply complementing the PDMC method 
with a re-sampling of the paths based on the value of the action, the re- 
sulting RQMC algorithm can afford system sizes typical of current BDMC 
simulations of continuous strongly interacting systems. Unbiased estimates, 
static responses and some insight into dynamical properties can be read- 
ily obtained. The dependence of the computational effort on the number of 
particles and on the quality of the trial function remains to be investigated. 
Such an analysis will determine whether this method can be useful in more 
general situations than ^He bulk liquid. Clusters, films and superfluids in 
restricted geometries are natural candidates for further applications. 
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For Fermion problems one has either to resort to the fixed-node approx- 
imation [32], or to cope with the sign problem [36]. In the former case, the 
dynamical information contained in the path is incorrect [4], but still the 
explicit knowledge of the weights along the path makes the algorithm free 
from the mixed distribution and the population control biases; furthermore 
it gives easily access to interesting quantities, such as a low-variance esti- 
mator of Born-Oppenheimer forces in electronic systems [37]. In the latter 
case the dynamics is correct, and can be used for example to get informa- 
tion on the ground and excited states from the imaginary-time evolution in 
the transient regime [38]; in similar cases however the real bottleneck will 
remain the sign problem. 

RQMC, VPI [27] and the technique discussed in Ref. [30] sample ex- 
plicit expressions of the imaginary-time propagator with the Metropolis 
algorithm to calculate zero temperature properties. We believe that these 
methods are very promising and deserve more attention than they have re- 
ceived so far. All these methods are based on the Metropolis algorithm: they 
enjoy therefore of a large freedom in the choice of the transition probability 
[39], which is probably not yet fully exploited. 
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